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ABSTRACT 

The problem of a porous graphite fiber plate subject to 
combustion is formulated. The transient one-dimensional 
model leads to two heat transfer equations on the graphite 
and air, and a mass transfer equation on the oxygen. In 
addition to a combustion heat generation term, the heat 
transfer model includes the mechanisms of conduction, radia- 
tion between fibers, and convection due to induced air flow 
through the mat. The mass transfer model includes diffusion 
and convection mechanisms, as well as a combustion consumption 
term. The temperature dependency of the system parameters 
is taken into account. The resulting nonlinear, coupled 
partial differential equations are solved by a Galerkin formu- 
lation of the finite element method. A number of computer 
analysis results are presented. The results show the effects 
of initial conditions, plate thickness, fiber diameter and 
air flow rate on ignition and extinction. 
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I. 



INTRODUCTION 



The combustion of a composite laminate plate consisting 
of graphite fibers in an epoxy matrix is discussed. The 
combustion process for a graphite laminate takes place in 
two stages, epoxy combustion followed by combustion of the 
graphite. At the end of the first stage the epoxy has burned 
away exposing a porous graphite mat. Spacing between the 
fibers is maintained by the residue of the combustion pro- 
ducts. This present work is concerned with the last stage 
of the combustion process, specifically, the thermal response 
of the remaining graphite fiber mat. 

The selection of the second stage for the analysis was 
made as a result of observations during composite plate 
"burn" experiments at the Naval Weapons Center, China Lake [1] . 
The objective of the experiments was to assess damage of com- 
posite aircraft structures subjected to open deck fires on 
aircraft carriers. For the experiments, an epoxy-graphite 
laminate was mounted into the wall of a wind tunnel parallel 
to the flow. One side of the composite plate was exposed to 
the wind tunnel flow and the opposite side was open to the 
environment. This arrangement was to simulate a composite 
wing section or fuselage subjected to typical wind condi- 
tions present on a carrier flight deck. It was observed 
that after being subjected to a fire which resulted in the 
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epoxy burning away, the graphite fibers would continue to 
smolder or burn. In other instances, the fibers would cool 
depending on the flow velocity and plate thickness. To 
understand this behavior better, a mathematical model was 
formulated to predict the conditions for which the exo- 
thermic process is self-sustaining. 

The air flow over one surface produces a pressure 
differential across the plate which induces a convective 
air flow through the porous medium governed by Darcy's Law. 

As a result, there is an enhancement of internal convection 
heat transfer, as well as a source of oxygen for combustion. 
The heat transfer mechanisms included in the model are 
(1) conduction, (2) convection, and (3) radiation. In addi- 
tion, non-volatile combustion is included in the energy 
balance as a heat generation term of Arrhenius type. The 
oxygen concentration is accounted for by a molecule-mass 
balance which includes (1) molecular diffusion, and (2) mass 
transfer by convection. Consumption of oxygen due to combus- 
tion is accounted for by a term similar to the heat generation 
term. 

The formulation of a one-dimensional model is presented. 
All properties are treated as temperature dependent in the 
transient analysis. To account for heat transfer mechanisms 
between the air and porous graphite medium, the air and 
graphite temperatures are treated as independent variables. 

The energy balance on the graphite and on the air, and the 
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molecule-mass balance on oxygen result in a system of three 
nonlinear, coupled partial differential equations. 

Integration of the field equations is accomplished by 
a Galerkin formulation of the finite element method. Due 
to the inherent stiffness of the field equations in the time 
domain, a modified implicit-Gear integration scheme developed 
by Franke [2] is used for a more efficient solution. 
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II. THEORY AND BACKGROUND 



A. INTERACTION OF HEAT TRANSFER AND COMBUSTION 

In this work, the combustion model proposed by Semenov 
was adopted. It is described in the texts of Frank-Kamenetskii 
[3] and Vulis [4] . A brief discussion of those features of 
the model which relate to the present investigation will be 
given. Fundamental to the model is the relation of reaction 
rate to temperature, and the interaction of the heat genera- 
tion and the heat transfer from the system. In general, the 
reaction rate, r, may be shown by the Arrhenius law as, 

r = A(T,tf>) exp(-E/RT) 

where A ^ is the characteristic time of the chemical reac- 
tion, E is the activation energy, R is the universal gas 
constant, T is absolute temperature and $ is the concentra- 
tion of oxygen. The concentration appears as cf> n in A ( T , ) 
for an nth order reaction. The heat generated by the exo- 
thermic process is obtained by multiplying the reaction 
rate, r, by the enthalpy of formation of the combustion 
products . 

In Figure 1, the heat generation, q , is plotted as a 
function of temperature. The curve is referred to as the 
S-curve for apparent reasons. The S-curves are distinguished 
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by two regions. Lower temperatures and reaction rates 
characterize region I, while region II is characterized by 
higher temperatures and reaction rates. In region I, for 
low reaction rates, there is an abundant supply of oxygen. 

The reaction is controlled by the temperature and the region 
is known as the kinetic regime. Here, the reaction rate 
increases exponentially with increasing temperature. In 
region II, the reaction is limited by oxygen concentration 
and is known as the diffusion regime. The reaction's weak 
dependence on temperature in region II is shown by the char- 
acteristic flattening of the S-curve in Figure 1. For higher 
concentrations of oxygen, the diffusion region is associated 
with higher reaction rates. The interaction of heat genera- 
tion, q^, and heat loss, q^, will now be described. 

In addition to the S-curve, Figure 1 shows two heat 
loss lines, q.^ and q They represent the heat transfer 

by convection for air temperatures, T^ and T 2 , respectively. 

The intersection of the heat generation curve and the heat 
loss curve, q^ = q^, defines the stable quasi-stationary point, 
A. As the air flow temperature increases from T^ to T 2 , the 
quasi-stationary temperature increases continuously from T^ 
to Tj. For temperature T ^ , q^ is tangent to q^ and an 
infinitesimal increase in the air flow temperature results 
in a jump of the reaction temperature from T_ to T_. At point 
I, which is defined as the "critical ignition condition", 
the reaction moves from the kinetic to the diffusion regime. 
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Temperature Tj is referred to as the ignition temperature. 
Frank-Kamenetskii states the ignition temperature of graphite 
is approximately 1100 to 1300 degrees Celsius. In order to 
keep the reaction model simple, the present investigation was 
limited to the kinetic regime. Although the description of 
the heat transfer process just presented is for convection 
only, the underlying ideas are valid for the conduction and 
radiation heat transfer mechanisms as well. The slope of the 
heat loss lines, , is equal to the product of the internal 
heat transfer coefficient and the surface area of fiber per 
unit volume. The slopes would not appear constant, and slight 
curvature in the heat loss lines, q , is apparent when radia- 
tion effects are included. As shown in Figure 1, an increase 
in the heat transfer rate resulting, for example, by an 
increase in the internal flow rate would result in a lower 
graphite temperature for a given air temperature. 

One reason for limiting the model to the kinetic regime 
was due to the jump in temperature and reaction rate at 
ignition. Though the discontinuity may be accounted for, the 
results show that for the region of temperatures approaching 
the diffusion regime, the heat generation will dominate the 
process. The system temperatures then rise rapidly. Other 
factors that limit the model to the lower temperatures are 
(1) the chemical reaction fixing the heat of combustion (com- 
plex at higher temperatures) and (2) the rapid consumption of 
fibers after ignition. 
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B. HEAT TRANSFER EQUATIONS FOR POROUS MEDIA 

In their investigation of oil recovery from underground 
reservoirs. Green and Perry [5] developed heat transfer 
equations for the solid and fluid phases in a porous medium. 
Their model included the three basic mechanisms: (1) physical 

movement of the fluid which carries its own heat capacity, 

(2) conduction of heat through the solid and fluid phases, 
and (3) convective heat transfer between the solid and fluid 
phases. Radiation heat transfer was assumed to be negligible, 
and combustion was not a consideration. The differential 
equations for the fluid and solid phases 

2 

3T f 3T f 3 T f 

Fluid: Pf c fP~ 9 t~ = “ Up f C f p Tx~ + k f P 2 hz(T f -T g ) (II. 1) 

3 x 

3T 3 2 T 

Solid: p s c s (l-p) T |- = k s (l-p) f+hz(T f -T s ) (II. 2) 

3 X 

were solved numerically by the finite difference method. 

In equations II. 1 and II. 2, c is the specific heat at constant 
pressure, u is the fluid velocity through the medium, h is 
the coefficient of heat transfer, and k is a pseudothermal 
conductivity. The model did not include change of properties 
with temperature. Their numerical results showed good agree- 
ment with their experimental results. Two useful conclusions 
derived from their analysis were: (1) both conduction and 

convection heat transfer mechanisms are important for internal 
Reynolds numbers less than one, and (2) for values of the 
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1/2 

dimensionless parameter, £ = (hz/k^) ' (k^/p f c^u) greater 
than .342, the fluid temperature is approximately the solid 
temperature. Riaz [6] proposed the following equation 
resulting from the equal fluid-solid temperature assumption. 



P 



3T , 
S C s3t 



3T 

P f C f U 3^ 




(II. 3) 



where the coefficients are defined as those in equations 
II. 1 and II. 2. For the present model, the more general 
approach ( i .e . , Green and Perry) will be taken in the formula- 
tion of the energy equations . 
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III. FORMULATION OF THE FIELD EQUATIONS 



A. INTERNAL FLOW MODEL 

Referring to Figure 2, the porous mat may be represented 
by a flat plate. One side is exposed to still air, and the 
other to an air flow of velocity, U ot . The environment on 
both sides is at ambient tempera ture, T^. 

The first step in developing the model was to derive 
an expression for the flow through the porous plate. Effects 
of the combustion process on the physical properties were 
neglected. Specifically, these were (1) variations in the 
porosity due to fiber consumption, and (2) effects on the 
viscosity and density of air from the introduction of combus- 
tion by-products into the flow. Preliminary calculations 
showed that the porosity and density of the material did not 
change appreciably until temperatures reached approximately 
2000 degrees Fahrenheit. Further, since eighty percent of 
air by volume is essentially inert, it was reasonable that 
its properties would not be significantly altered during the 
combustion process. Asstiming a heterogeneous material of 
constant porosity, a Reynold's number for the interior flow 
is defined as 



Ee i ■ P a u p d/,J 



(III.l) 
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where is the mass density of air, is the pore velocity, 
d is the diameter of the fiber, and u is the dynamic viscosity 
Extensive experimental work shows that Darcy's law for flow 
through a porous medium is valid for a limited range of 
Reynolds numbers. Muskat [7] proposed that Darcy's flow is 
valid for Re^ < 1. Scheidegger [8] points out that a number 
of investigators give a much higher upper limit. Preliminary 
calculations of Re^ for the model showed typical magnitudes 
on the order of 1. Neglecting the gravity term, Darcy's law 
for flow in porous media is. 



u 

P 



m dP 
y dx 



(III. 2) 



assuming a constant pressure gradient, equation III. 2 
becomes , 



u 



P 



m A£ 
y L 



(III. 3) 



where AP is the pressure differential across the plate 
thickness, m is the specific permeability of the medium, 
and y is the dynamic viscosity. In equations III.l and III. 3, 
p and y were treated as temperature dependent. 

Specific permeability is a measure of a porous medium's 
ability to allow fluid to flow through. It is dependent only 
on geometry and the physical nature of the medium. It has 
dimensions of length squared. Specific permeability is often 
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called the hydraulic conductivity due to the similarity of 
Darcy's law to that of Fourier's law for heat conduction. 

The permeability for a particular medium is normally measured 
by experimentation. However, there are several empirical 
models that may be used to obtain values for permeability that 
are in agreement with experimental results. For the idealized 
geometry of the porous medium shown in Figure 2, the permea- 
bility was obtained by a serial type model, equation III. 4, 
proposed by Scheidegger [8] , 



and D is the thickness of a ply. Porosity has units of 
volume of void per unit volume of medium. The average pore 
diameter, <5 was defined as. 



for the geometry in Figure 2. The tortuosity, t, is a measure 
of length of travel for a fluid particle per unit thickness 
of the medium. Referring to the geometry of Figure 2, the 



m 




(III. 4) 



where the porosity, p, is defined for a prous material 



comprised of fibers with cylindrical shape as. 



P 



1 - |-(d/D) 2 



(III. 5) 



6 



(s + D)/2 
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tortuosity depends on the ratio, d/D. The lower limit 
occurs when d/D equals zero and for the upper limit, d/D 
equals 1. Carman [9] presents a tab.le of tortuosities for 
given geometries and particle shapes. For the model. 

Carman's suggested value of t = 1.4 was used. 

To obtain the pressure differential, AP, across the 
porous plate, Bernoulli's equation was used. The following 
observations showed this to be a valid assumption. Schlichting 
[10] points out that for the ratio of u^/U^ in the range of 
.0001 to .01, the effects of "blowing" or "suction" on the 
potential flow over the porous plate may be neglected. A 
typical value of u /U for the model at which = 25 kt. was 
.0028. For steady flow over a flat plate, the flow field 
outside the boundary layer may be described by Bernoulli's 
equation. This is a direct result of the Navier-Stokes equa- 
tion. In addition, the pressure gradient across the boundary 
layer may be taken equal to zero. Therefore, Bernoulli's 
equation (eq. III. 6) may be used to obtain the pressure 
differential across the plate, 

2 

, p u 

I — dx + -s- = constant. (III. 6) 

J p a 2 

The parameters in equation III. 6 are defined as follows: 

P, pressure; p , the density of air; U, velocity of the air 

cl 

over the flat plate. For the model, the density of air was 
approximated by the Ideal Gas law as 



23 



(III. 7) 



p = P/R T 
a a a 



where R is the gas constant for air, T is the absolute 

3 . 3 

temperature of air, and P is the pressure. Substituting 
this into equation III. 6, gives 



R T 



/ ,-JLJL,dP + 



= constant 



and upon integrating, the expression becomes 



R T In P + Hr- 
a a 2 



= constant 



or, 



U 1 

R a T al ln P 1 + ~ 



°2 

R a T a2 ln P 2 + 2 



Letting 



P, = P , U, = 0, T = T _ = T , P=P. U, = U 
1 1 al a2 00 2 l 2 <> 



and siobstituting these in above, yields 



u : 

R T In P = R T In P T + - T 

a » 00 a 00 L 2 



Rearranging the previous expression gives, 
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R T In (P /p_ ) = -£■ 

a » « L 2 

Taking the exponential of both sides and solving for P T 

JL 

results in 



P_ = P exp (-U 2 /2R T ) 

L 00 00 a 00 

From the above expression and noting that AP = P^ - P^, 

AP may be expressed as 

AP = P [exp (-U 2 /2R T ) - 1] 

00 00 a 00 

Darcy's law in approximate form (eq. III. 3) for the model 
is 



u 

P 



m AP 
V L 



and substituting for AP, it follows that. 



u 



<” ^[l-exp(-uf/2R a ij]) 



(III. 8) 



where 



m 






25 



The pore velocity, u of equation III. 8 is for an 

poo 

isothermal medium at ambient temperature, T . Since the 
model is to investigate the transient problem, a general 
expression to obtain pore velocity at different temperatures 
had to be developed. 

Consideration of the continuity equation for one-dimen- 
sional steady flow. 



d (p u )/dx = 0 

a p 



(III. 9) 



becomes 



u (dp /dx) + p (du /dx) = 0 

p a a p 



or 



u (dp /dx) = - p (du /dx) 

a. d 

Multiplying through both sides by dx, the expression 
becomes , 



u dp 
P 



a 



p 



a 




Separating variables and integrating both sides, gives 
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J 



p 



CO 



da) 

to 



u 




u 

poo 



dco 

to 



and the integral becomes 



ln(p /p ) = ln(u /u ) 

a °° p OT/ p 



Taking the exponential of both sides, gives 



p./p b = u /u 
a 00 p°° p 



Rearranging, the pore velocity is expressed as 



u = p u /p (III 

p p« a 

Referring to the Ideal Gas law for the density of air, and 
substituting this into equation III. 10 gives 



u 



= (p 



u R T )/P 
= p°° a a 



(III .11) 



The pore velocity is now a function of air temperature and 
pressure. However, in preliminary calculations, the pressure 
difference, AP was small and for the range of parameters of 
interest, it did not exceed five percent of the magnitude of 
ambient pressure. Therefore, as a simplification, it was 
assumed that P equals the ambient pressure in equation III. 11, 
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and the expression becomes 



u 




Noting that p ra = P^/R T^, the pore velocity may be expressed 
as 



where u is obtained from equation III. 8. T and T are 

poo ^ a 00 

the absolute temperatures of the air at a point in the interior 
of the medium and of the environment, respectively. 

B. ENERGY EQUATIONS 

Here we consider the development of the energy equations. 
However, before continuing, discussion must be made as to 
the approach taken for energy balances in porous media. In 
previous work by Denbigh and Turner [11] , and by Co 11a day 
and Stepka [12] , it was suggested that the temperature of 
the porous solid and of the fluid be considered equal. This 
greatly simplifies the formulation of the problem. However, 
under certain conditions, as shown by Green and Perry [5], 
the assumption of equal temperatures may not be valid. With 
recent developments in numerical techniques permitting the 
solution of non-linear systems, the development of the present 
model for heat transfer in a porous medium was not restricted 



u 



P 



u T /T 

poo a ' 00 



(III. 12) 
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to the assumption of equal temperatures. Therefore, the 
temperatures of the air and of the porous solid are treated 
as independent variables in the present formulation. 

The model is based on two assumptions. First, due to 
the small thickness of the graphite fibers, the temperature 
across each individual fiber was assumed to be constant. 
Second, for a similar size consideration, the temperature of 
the air in each individual pore was assumed constant. How- 
ever, the temperature from fiber to fiber and from pore to 
pore was not restricted, and was allowed to vary according 
to the heat transfer mechanisms described next. 

To perform energy balances on both the porous solid and 
on the air, a differential volume of porous material was 
segregated into respective volumes of constituents, that is, 
graphite fibers and air (Figure 3) . An energy balance may 
be done on each volume separately. The convention used for 
the balance is: 



Heat into + Heat 

dV Generation 



Heat out 
of dV 



Increase in 
Internal 
Energy 



where 



dV = dx dy dz 
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1 . Fiber Heat Transfer Equation 

Considering the differential volume, (l-p)dV of 

graphite fibers, and "smearing" the fibers to form a macro- 

scopically homogeneous material, the heat transfer mechanisms 

are shown for the x-direction in Figure 4. The heat transfer 

mechanisms are: q conc j' heat conduction through the medium; 

q , convection from the fibers to the air; q . , radiation 
^conv ^rad 

transfer from fiber to fiber; AHr(Tg,<j>) heat generation rate 
per unit volume. AH is the enthalpy of formation, and 
r(T ,<{>) is the fiber mass consumption rate. This term will 
be discussed in detail later. 

Representing the terms in Figure 4 by Taylor series 
expansions and combining them with respect to the convention 
previously stated, the energy balance is 



q cond + q rad + 4Hr <V*> (l-P)dxdydz = q cond 



+ dx 
3x 



3 q . 

+ q ■, + dx + q + -^(l-p) dxdydz 

^rad 3x M conv 3t * J 



Subtracting terms, q cond and q rad from each side and 
rearranging, the expression becomes (neglecting the higher 
order terms) 



3q 



cond 



3x 



dx - 



aq 



rad 



3x 



- q + AHr (T , ) (1-p) dxdydz 

^conv g J 



3 E 
3 1 



(1-p) dxdydz 



(III. 13) 
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The heat transfer terms in equation III. 13 will now be 
considered. 

a. Thermal Conduction 

Fourier's law for conduction may be written for 
the incremental volume as 



q , = - k ^dydz 

^cond g 3x ■* 



(III. 14) 



where k^ is the effective conductivity of the porous solid, 
and Tg is the temperature of the fibers. There are a variety 
of models for the effective conductivities of porous media. 
They depend largely on the geometry and on the nature of 
the constituents. For the cylindrical shape of the fibers, 
the empirical relation proposed by Rohsenow [13] was used 



k = k {1-&[1 + 



Y i^/l- (1~V ) 



2 , 1/2 



(1-Y) 



1/2 



ln(- 



-)]} (III. 15) 



2 

for y <1. The parameters of equation III. 15 are defined 
as follows: 3 = d/D (refer to Figure 2) , y is 



= it. 



A 



(1-A) 



] 



where A = k /k', k is the conductivity of the air; and k' 
a g a g 

is the actual conductivity of the graphite in bulk form. 
Since k^ depends on temperature, k^ was also treated as 
temperature dependent. 
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b. Radiation Heat Transfer in Fibers 



Radiation heat transfer in the model was repre- 
sented by an analog to Fourier's law of conduction as 



^rad 



3T 

3x 



dy dz 



(III. 16) 



where is an equivalent conductivity of the fibers due 
to radiation. Simple assumptions and approximations may be 
made to show this. Assuming air to be transparent to radi- 
ation, and treating a ply of graphite fibers as an infinite 
wall, the net heat flux between plies may be written as 



*rad 



— (T^ 

(2-e) gx 



— T ) 

gx+dx' 



(III. 17) 



where and T gX+( j x are the respective absolute ply tempera- 

tures, e is the emissivity of the graphite, and a is the 

Stefan-Boltzman constant. Treating T as the variable 

gx 

temperature, q” a ^ may be expanded in a Taylor series about 
T gx+dx' Ne 9 lectin 9 the higher order terms, the series 
expansion may be written as 



*rad 



oe -(T 4 



- i 4 . . ) + i 3 



(T - T 



2-e ' gx+dx gx+dx 2-e gx+dx gx gx+dx 



) 



Simplifying, the above equation becomes, 



4 rad 



— — - T^ (T - T ) 

2-e gx+dx gx gx+dx 



(III. 18) 
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Taking Fourier's law for heat transfer, 



q 



II 

rad 



k 

r 



dT 




dT 



and approximating with 





T 



gx+dx 

Ax 



A 




q 



it 

rad 



becomes 



q 



ti 

rad 



k 



r 



(T , . - T ) 

gx+dx gx 

Ax 



(III. 19) 



Multiplying the right side of equation III. 18 by Ax/Ax, 
equation III. 18 and equation III. 19 may be equated 

/V /V /s /s 

(T - t ) (t — - T ) 

'• . gx+dx gx _ _ 4a eAx ,13 gx+dx gx 

^rad r Ax 2-z gx+dx Ax 



Cancelling terms, k^ becomes 

k = 4aeAx ^,3 
r 2-z gx+dx 



(III. 20) 



where ax is now equal to <$ , the average pore diameter. 

From the close spacing of the fiber layers, the temperature 
difference will be small as compared to the magnitude of the 
temperature. Noting this the average absolute temperature of 
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The equivalent 



the fibers, T may be substituted for T , , . 

g J gx+dx 

radiation conductivity expression becomes 



k = ■ T 3 (III. 21) 

r 2-e g 

and the radiation heat transfer from fiber to fiber may be 
represented by, 



^rad 



k 



r 




4ae6 ,13 
2-e l g 



3T 



3x 



dy dz 



(III. 22) 



If the temperature gradient in the fibers is not 

large, then equation III. 22 will be a good approximation of 

the radiation heat transfer between the fibers. 

c. Convection Heat Transfer 

The convection heat transfer term in Figure 4 

was treated in a similar manner to that of the one-dimensional 

fin equation. q was introduced as 

conv 



q 



conv 



h. (T - T ) dA 
l g a 



(III .23) 



where T^ is the temperature of the graphite, T^ is the 
temperature of the air, dA is the surface area of graphite 
in the differential volume, and h^ is the internal convection 
heat transfer coefficient. An empirical expression in the 
form of a Colburn j -factor was developed by Yoshida, 
Ramaswami, Hougen [14] . This was used to determine h^. 
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h • C )i T / *3 Cl 

J = — -i-[^_] z/J = .91 Re * ”* t q (III. 24) 

c b o a 

for Re < 50. is the effective specific heat of the porous 

medium and is defined as, 

c b = pc a + (1 ' p)c g 

where c^ is the specific heat of air at constant pressure; 

c^ is the specific heat of graphite; and p is the porosity 

as defined by equation III. 5. G q is a pseudo mass velocity 

defined as G = p u p. Also in equation III. 24, y is the 
o a p r 

viscosity of air, and k is the conductivity of air. The 

3 . 

Reynolds number appearing in equation III. 24 is defined as 

Re = G /z y f 
o o 

and is not the same as Re^ defined previously for Darcy's 
law. z is the surface area of graphite fibers per unit 
volume, and is determined from geometrical considerations 
which will be described shortly. The parameter, , is a 
dimensionless shape factor which depends on the geometry of 
the fibers. For a cylindrical fiber shape, V is equal to 
.91. From the assumption of constant porosity, z remains 
constant. However, the air properties are temperature dependent; 
therefore, h^ is allowed to vary with temperature. 
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d. Heat Generation Rate 



The reaction rate, r(T , <J>), whether based upon 

g 

the Collision Theory or on the Theory of Absolute Reaction 
Rates, results in an expression of Arrhenius type, 

r = A (T , <j> ) exp (-E/RT) 



where E is the activation energy, R is the universal gas 
constant, and T is absolute temperature. Although E is well 
defined for the combustion of graphite, A(T,<j>) is a function 
of temperature and concentration and depends upon chemical 
kinetic theory. Development of a general model for the chemi- 
cal kinetics is beyond the scope of this work. A less ele- 
gant, but useful approach of using a relation obtained 
experimentally was taken. A relation derived from the work 
of Parker and Hottel [15] was used to predict the combustion 
rate of graphite fibers. 



= 4 . 014 x 10 R, T 



1/2 



y <P exp ( 



-39,883 



) 



Ibm 

ft 3 -hr 



(III. 25) 



The development of equation III. 25 from the original work 
may be found in the appendix. This expression assumes a 
simple first-order reaction for the combustion of graphite 
in air. However, Frank-Kamenetskii [3] has further refined 
Parker and Hottel' s expression to yield for the present 
model. 
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r 



lbm 



(III. 26) 



4.014 x 10 




-39,883 



T 

g 




The combustion is now characterized as an nth order react- 
tion, where n is in the range of 1/3 to 2/3. This form of 
the reaction rate expression better approximates Parker and 
Hottel's experimental data. For the analysis, n was chosen 
to be 1/2. However, it should be pointed out that the behavior 
of the reaction varies significantly within the range of n 
proposed by Frank-Kamenetskiii . 



combustion, it was assumed that the dominant chemical reac- 
tion is 



This is an idealized treatment of the reaction that actually 
takes place. Kanury (16] states that, the actual reaction 
is more complex, yielding various concentrations of carbon 
dioxide and carbon monoxide, depending on the temperature. 
Noting the abundant oxygen supply present in the kinetic 
regime, it is reasonable that the leaner oxide is the preva- 
lent by-product of combustion. The fuel-oxygen ratio, f, 
for this mechanism is 12/32. Carbon monoxide may be formed 
when there is little free-oxygen present, such as, during 
combustion in the diffusion regime. The stoichiometric 



In limiting the model to the kinetic regime of 



C + 0 



2 



CO 



2 
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fuel-oxygen ratio, f = 12/16, for this mechanism is naturally 
greater. To obtain the heat generation rate, the reaction 
rate, r, was multiplied by the enthalpy of formation, AH, 
of carbon dioxide. AH has units of Btu per lbm of graphite 
consumed. 

e. Change in Internal Energy 

The last term of the graphite energy equation, 

d E 

the rate of change of internal energy, -r-r-, results by 

O t 

setting E = p c T . Thus, 

g g g 

3 T 

l^-d-p) dxdydz = (1-p) P g c g dxdydz (III. 27) 



where p^ is the density of graphite in bulk; c^ is the 
specific heat of graphite; and p is the porosity. 

The heat transfer mechanisms for the graphite 
fibers, equations III. 14, III. 16, III. 23, III. 26, III. 27 are 
substituted into equation III. 13 which upon dividing through 
by dV = dxdydz yields the heat transfer equation for the 
fibers , 



( 1 -p) f-(k 



Ti? 5 + ( 1 " p) k (k r 


9T 

— SL) 

9X 


- h. §^(T - 1 
1 dV g 


+ (l-p)AHr(T , 4 .) 

y 


= 


9T 

(1 -P> p g C 9 Tf 



The coefficient, (1-p), accounts for the fact that not all 
of the volume of porous medium is comprised of graphite. 
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Dividing through by (1-p) , combining the first two terms, 
and defining dA/dV as z, the surface area of graphite per 
unit volume, the expression becomes 




+ k 



3 T 

) — 2 ] 
•' 3x J 



TT^T'Tg - V + 4 Hr <V+) 



p 



c 

g g 




(III. 28) 



The parameter, z, is obtained from geometry considerations. 
For a cylindrical fiber shape, the expression for z is. 



z = -rd/2D^ 



(III. 29) 



The one-half factor in equation III. 29 accounts for the un- 
certainty of the actual amount of exposed fiber surface area, 
which resulted from the initial combustion of epoxy. The 
non-dimensionalization of equation III. 28 is presented in 
Appendix A. 

2. Internal Flow Heat Transfer Equation 

The energy balance on the internal flow is obtained 
by a procedure similar to the energy balance on graphite. 
Figure 5 presents a "smeared" differential volume of air, 
pdV, which shows the heat transfer mechanisms. The heat 
transfer mechanisms are: 9 con( j' heat conduction through the 

A 

air; q , convection from the fibers to the air; m h, 
conv ci 
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energy transport due to the air flow, where h is the 
enthalpy of the air . 

Performing the energy balance (neglecting the higher 
order terms) yields. 



*cond + m a k + q conv q cond 



3q 



cond 

3X 



dx 



• , 3h , 3E , , , 

+ mh + m — dx + TT- pdxdydz 
a a o x o c 



(III. 30) 



Cancelling terms and rearranging, the expression becomes, 



3q cond j • 3h j 
— — dx - m a 3^ dx + q conv 



pdxdydz 



(III. 31) 



a. Thermal Conduction 

As before, q conc j was replaced by Fourier's law 
of heat conduction in the form, 



3T 

q cond ~ k a 3x d ^ dz 



(III. 32) 



where is the conductivity of air. 

b. Energy Transport by Flow 

Using the perfect gas assumption for air, the 
enthalpy h was expressed as. 



dh = c dT‘ 
a a 
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where is the specific heat of air at constant pressure. 
The specific heat is treated as constant since it does not 
vary appreciably over the range of temperatures associated 
with the problem under investigation. It follows that, 



ah 

ax 



3T 

a 

c a ax 



(III. 33) 



Assuming the effects of the combustion by-products on air 
are negligible, continuity considerations give the mass 
flow as. 



iti = p u^ dydz (III. 34) 

0 0. P 

oo ^ oo 

where the density of air and the pore velocity are evaluated 

at ambient conditions . 

c. Convection Heat Transfer 

The convection heat transfer, q given by 

n conv 3 1 

equation III. 23, with h^ obtained from equation III. 24, 
becomes 



^conv - h i ( V T a )dA (HI. 35) 

d. Change in Internal Energy 

For a perfect gas, the internal energy is, 

dE = p c dT 
a a a 
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where is the specific heat of air at constant pressure. 

From this it follows that 

9 T 

p dxdydz = pc p dxdydz (III. 36) 

d U a a oX 

Substituting equations III . 32-III . 36 into equation 
III. 31, the heat transfer equation for the internal flow 
becomes. 




3T 

•j — ) pdxdydz 



3T 

™a c a Tx P dxd y dz 



3T 

+ h • dA (T - T ) = pc p dxdydz 

l g a a 3t c 



As in the energy equation for the fibers, the 
porosity, p, accounts for the fact that some of the volume of 
the porous medium is not air. Dividing the last equation 
through by pdxdydz and letting z denote dA/ dxdydz, yields 




3 T 

a 



x 



m c 
a p 



3 T 
a 

3x 



h . z 
i 



(T - T ) 
g a 



P 



c 

a a 



3T 

a 
3 1 



(III. 37) 



Non-dimensionalization of equation III. 37 may be found in 
Appendix A. 
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3 . Oxygen Transport Equation 



The final consideration in the formulation of the 
model is the mass transport of oxygen. The mass transport 
equation is obtained by a mass balance on a differential 
volume of porous medium. In accordance with convention, 
the mass balance is given by 

Oxygen into _ Oxygen out + Oxygen Oxygen 

dV of dV Consumption Accumulation 



Figure 6 presents a differential volume, pdV, showing the 
relevant molecule-mass transport mechanisms. The mass 
transport mechanisms are: mg is the molecular diffusion; 

m c is the convective transport due to mass flow; ^r(T^,<}>) 
is the consumption of oxygen due to combustion; r is the 
reaction rate, and f is the stoichiometric ratio of the 
reaction. Representing the terms in Figure 6 by Taylor 
series expansions, the molecule-mass balance becomes (neg- 
lecting higher order terms) , 



ra B + m c 



3m B 

"B + 3X dX + + 



8m 

c 
8 x 



dx 



+ ir ( T , <f> ) pdxdy dz + f-^-pdxdy dz 

I CT d tl 



(III. 38) 



Upon cancelling terms and rearranging, equation III. 38 
becomes 
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II- pdxdydz (III. 39) 

a t 



3m 3m . 

TF dx _ TF dx " I r (T g' <t ‘ ) P dxd Ydz = 



where 4> is the concentration of oxygen. The individual mass 
transfer terms in equation III. 39 will now be developed, 
a. Molecular Diffusion 

Molecular diffusion m^, is given by Fick's law, 



m 



B 



- B |-^ dy dz 

3 X 1 



(III. 40) 



where B is the diffusivity of oxygen into the porous medium. 

The diffusion of gases in porous media is limited 
by two factors. One is the collision of the gas molecules 
with other molecules, and the second is the collision of the 
gas molecules with the solid pore walls. This second phenom- 
enon is discussed by Bennett and Myers [17] and is referred 
to as Knudsen diffusion. To establish which behavior dom- 
inates, the mean free path of the gas (in this case, oxygen), 
must be calculated. Using the expression given by Treybal 
[18] , 



w = (3.2y/P) (RT/2irg M) 1/2 

where g is the gravitational constant, the value obtained 
c 

for mean free path, w, is smaller than 0.16. Thus, the 
collision of the oxygen molecules with those of air will 
restrict the diffusion process. On the contrary, if the pore 
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diameter is smaller than the mean free path, then the colli- 
sion of the molecules with the solid walls will be the 
limiting factor. In the intermediate case, the mean free 
path and the pore diameter are of the same order of magni- 
tude. In this case, Bennett and Myers [17] proposed the 
following expression be used to obtain a value of diffusivity. 




(III. 41) 



where B R is the resulting diffusivity; B^ is the diffusivity 
based on molecule-molecule collisions; and B^ is the diffusivity 
based on molecule-pore wall collisions . 

For the problem under investigation here, prelimi- 
nary calculations of the mean free path of the oxygen showed 
that diffusion is primarily of the inter-molecular collision 
type. With respect to this, and knowing that molecular 
diffusion is highly temperature dependent, the diffusivity, 

B* was obtained from an expression presented by Gilliland 
[19] . 

m3/2 

B* = 435.7 (V^ /3 + vJ/ 3 )" 2 (M" 1 + M~|) 1/2 cm 2 /s 

(III. 42) 

where P is the absolute pressure; and V ^2 are the molecular 

volumes of air and oxygen respectively; M and are the 

3 U 2 . 

molecular weights of air and oxygen, respectively, and B* 
is the diffusivity. For this expression, T is in degrees 

3 
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2 

Kelvin, and B* has the units of cm /s. As a result of 
tortuosity and porosity affecting the process, Denbigh and 
Turner [11] proposed that the effective molecular diffusivity 
for a porous medium is given by 



B e = pB*/ t . 



The porosity p accounts for the void in a differential cross- 
section . 



by 



b. Convective Transport 

The convective term in equation III. 39 is given 



m c = pu p (f>dydz 



(III .43) 



The pore velocity u is a temperature dependent variable, 

P 

and is so treated. 

c. Consumption Rate 

The reaction rate term, r(T , <j>) was discussed 

9 " 

previously as equation III. 26. The consumption rate of 
oxygen is determined by multiplying r by 1/f. 

Substituting equations III. 40 and III. 43 into 
equation III. 39, the oxygen transport equation becomes. 



fe ( T- f£> Pdxdydz - fj(u p +> pdxdydz 



(■£■) r (T , <j>) pdxdydz 



pdxdydz 
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Dividing both sides by pdxdydz, and letting B = B*/t, the 
expression becomes 



fe< B !£> - - ( I )r(T g 



_ 3 4 > 

'*> ~ Jt 



(III. 44) 



Expanding the second term of equation II I. 44, 




P 



4 >) 



u 

P 





4 > 



and noting that u^ depends on from equation III. 12, the 
expression yields 




P 



40 



u 



P 




au aT 

3T ax 

a 



<P . 



(III. 45) 



3u / 3T is 
P a 



the oxygen 



obtained from equation III. 12. 

Substituting equation III. 4 5 into equation III. 44, 
transport equation becomes. 



i_( B 11 ) _ 

ax lB a x 



u 



a 4> 

ax 



au aT 

p a 

aT ax 



- (£)r(t ,#) 



= (III. 46) 



In contrast with the heat transfer equations, mass balance 
does not depend on porosity. However, porosity does influence 
the oxygen concentration by appearing in the boundary condi- 
tions. Non-dimensionalizat ion of equation III. 46 is pre- 
sented in Appendix A. 
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IV. BOUNDARY CONDITIONS 



It was intended that boundary conditions be imposed 
which simulate the composite plate "burn" experiments of 
Fontenot [1] . However, for a one-dimensional model, the 
boundary conditions can only be approximated at best. The 
difficulty is associated with the existence of momentum and 
thermal boundary layers on the exterior surface of the por- 
ous plate. Complications also arise from the "blowing and 
suction" effects which are caused by the flow through the 
plate . 

The following boundary conditions provide a reasonable 
model for the one-dimensional problem, 




h, (T - T ) + ae (T 4 - T 4 ) 

1 g « g oo 



at x = 0 (IV. 1) 




cre(T 4 - T 4 ) at x = L (IV. 2) 

g oo 



T 



a 



T 



oo 



at x = 0 



(IV. 3) 



3x 



a 



3T 

2 . 

3x 



at x = L 



(IV. 4) 



B = U p U " P< 0 at x = 0 



(IV. 5) 




at x = L 



( IV . 6 ) 
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Boundary conditions IV. 1 and IV. 2 account for convection 



and radiation heat transfer from the porous solid. With 

an influx of air at x = 0 , the heat transfer coefficient, 

h^, will depend on the magnitude of the flow entering the 

plate. Neglecting "blowing and suction" effects, several 

empirical expressions based on free convection were used to 

obtain a heat transfer coefficient at x = 0 . These expressions 

2 

yielded a range of values for h^ from 1 to 3 Btu/hr*ft *F. 

As an approximation, h.^ was taken as the average value, 2.0 
2 

Btu/hr»ft *F. However, choosing h^ as any value in the 
above range did not affect the results of the analysis. 

The magnitude of the forced convection heat transfer 
coefficient at x = L, h 2 , varies in a direction parallel to 
that of the external flow, U^. In addition, h 2 depends on 
the efflux of the gas at the surface. To simplify the analy- 
sis, h 2 was approximated by the relation for a smooth flat 
plate given by Holman [20] as 



for laminar flow, where Pr is the Prandtl number and L* is 
a reference length, arbitrarily taken as unity. For turbulent 
flow. 



h 



2 




(IV. 7) 



h 



2 



7 $ Pr 1/3 (.037 Re 



8 



850) . 



( IV . 8 ) 
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Kays [21] provides an alternate scheme for treating boundary 
layers and "blowing and suction" effects. 

Boundary condition IV. 3 accounts for the air entering 
the porous plate at ambient temperature. The Dankwerts 
boundary condition proposed by Riaz [6] which accounts for 
conduction-convection interaction, is a more reasonable 
simulation. 

The boundary condition at x = L for the air heat trans- 
fer equation appears to be a reasonable one. It follows the 
behavior observed in the interior of the plate and does not 
fix the temperature of the air at the wall. The boundary 
condition, 3T /3x = 0, used initially provided results simi- 
lar to those obtained using equation IV. 4. 

Boundary conditions IV. 5 and IV. 6 are Dankwert conditions 
[22] for flow through porous media.' A convincing discussion 
for the Dankwerts conditions is given by Bischoff [23] . A 
brief summary of the discussion is presented in Appendix D. 
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V. FINITE ELEMENT METHOD 



A. GALERKIN FORMULATION 

A Galerkin formulation of the Finite Element Method was 
used to obtain solutions of the field equations . A conven- 
ient form of the field equations was used in the formulation. 
These equations are presented in Appendix A. 
graphite energy equation 

d 0 

t? (x i If' - V s - r) + x 3 r( V*’ = x 4 “ (va) 

air energy equation 






v 3 (0 - D 



a_r 

at 



(V. 2) 



oxygen diffusion equation 



3 / 3$ \ 3< & 3 r * /m xX 

t — ( uj. t — ) - a). — — $-ijj.r(T ,<b) 

3n 1 3n 23n 33n 4 g r 



- 3 $ nr -3 \ 

= at (v - 3) 



where the coefficients. are defined in Appendix A. 

Field equations III. 28, III. 37 and III. 46 subject to 
boundary conditions, IV. 1 - IV. 6 , and initial conditions 
define the problem. The closed domain (0,L) was partitioned 
into (n-1) contiguous elements of variable length 5,^, i = 1, 
. . . , (n-1) . This defines an n nodal point model. The three 
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non-dimensional field variables 0, r, $ were approximated 
by 



e(n,t) = 4' 1 (n,t) = ^G i (n)0 ± (t) (V.4) 

r(n,t) = 4*2 (n ,t) = ^G i (n)r i (t) (V.5) 

*(n,t) = 4' 3 (n,t) = Ig l ( n)* ± (t) (V.6) 

where G^, for i=l, n, is a set of specified basis 

functions with local support, and the sets {0. , r . , 4>, ;i = 1, ... 

Ill 

n} are the solution coefficients to be determined. The 

were selected to satisfy the condition G^(nj) = 6 ^ j , where 

the kronecker delta 6. . is defined by 6. . = 1 for i = j, 

ID ID 

and = 0 for i ^ j. As a result, 0, r, and $ are the 

values 4^, 4 , 2 r V ^ at the nodal points (i.e., 0^(t) = 4'^(n^/t)) 
Linear interpolation functions (Figure 7) were used 
as the basis functions. These are the lowest polynomial 
functions which provide the necessary function continuity. 

As a measure of error, a residual function, r., is 

l 

defined for each field equation by 

r ± (n»t) = A ± (4') - 4' i i = 1, 2, 3 (V.7) 

where A^ denotes the spatial operator of the ith equation. 

For convenience, the following convention for differentiation 
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is adopted, 





^ - < >• 

9 2 ( ) _ , * ,, 

. 2 ' ’ 

9n 

3( } = (*) 

at ' ' 


For field 


equations V.l, V.2, and V.3 the residuals are: 
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n n 

U,(0,r) l G ! 0 • ] ' - X~(T) l G. (0. - T.) 
x ~ ~ i=l 1 x ~ i=l 1 1 1 

n 

+ * 3 r(0,$) - x 4 l G.Q^ CV.8) 

i=l 



r 2 * 


[ vi ( r) l G i r i 3 ' - v 2 l G|r ± + v 3 (r) l G i ( 0 i -r i ) 
x i=l 1 1 z i=l x i=l 1 1 1 




n 

- V.(r) l G.r. (V. 9 ) 

4 - ill 11 



II 


n n n 

I G! *.] ' - u 2 (r) l G!$ i -o 3 (r) l g.t.g.*. 

x ~ i=1 11 z i=1 11 J i=1 1111 




n 

- to 4 r(0,$) - Ur 7 G.<t>. (V.10) 

4 ' “ 5 i=l 1 1 
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where the coefficients multiplying the response variables 
are themselves functions of the response variables, and thus 
the equations are nonlinear. In accordance with the Galerkin 
method, the final system of ordinary differential equations 
was obtained by setting each residual, r ^ , orthogonal to 
each basis function, G^, that is 



1 

J G.r.dn = 0 i = 1,2, ...,n; (V.ll) 

0 1 3 . , - 

j — 1,2, ... ,n 

The 3n ordinary differential equations given by equations 
V.ll retain the character of the original set of partial 
differential equations. Thus linear field equations trans- 
form to matrix operators and nonlinear, coupled field operators 
become nonlinear, coupled algebraic operators. Incorporation 
of the boundary conditions resulted in 3n nonlinear, coupled 
ordinary differential equations 



# 

a ( t ) 4 , (e,r,$) + f ( t) = b ( t) i 



(V. 12) 



subject to initial conditions, where B is a 3n x 3n matrix, 
A is the operator associated with the field operator in 
expression V.7, and F(t) is an excitation vector. 

Adopting the convention. 



< G i> 






n 

l 

i=l 



G.'V . 

l l 
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and applying the operation of expression V.ll with an inte- 
gration by parts on the second order derivatives gives. 



{Gi> A 1 <G i >' {0} 



1 

0 



1 

X i / {G. } ' <G . > 'dn {0} 

1 0 1 3 



1 1 

- f {G . } <G . >dn (0 } + x 0 j {G. }<G .>d n (r} (V.13) 

* o 3 o 3 



1 1 

+ A 3 / {G.}dnr (0, $) = A 4 / {G . } <G . >dn { 0 } 

0 ~ 0 3 



{G i }v l <G j > ' 



1 

0 



1 

v. / {G.}'<G.>'dn{r} 

o 1 : 



1 1 

- V 9 / {G. }<G.> 'dn{r} + V- / {G. }<G . >dn {0} (V.14) 

* 0 1 3 0 x 3 



1 1 

— v_ j {G . } <G . >dri { r } = v, j {G . } <G . >dv { r } 

3 0 1 J 4 0 1 3 
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1 



{G. }u>, <G.> ' {$} 

i 1 3 



0 



1 

w, / {G . } ' <G . > ' dri { $ ) 

1 ; 0 1 3 




1 



/ 

0 



{G. } <G . > ' dn { 4* } 



1 

co-, / {G. }<G.>' {r}<G.>dn{$} 

J o 1 3 3 



1 

- a) . / {G . }dnr (0, $) 

4 0 1 




1 



/ 

0 



{G ± } <Gj >dn { $} 



(V . 15) 



The first term in each of the above expressions is a boundary 
term which permits the incorporation of natural boundary 
conditions which will be shown in Section V.B. Coefficients 
T, v, co are variable-dependent properties, and were taken 
as the average value of the properties over an element. In 
the limit, as the elements get smaller (i.e., n <*>) , the 
average values converge to the exact values. 

Inspection of expression V.13, V.14, V.15 shows the five 
operators , 



1 

/ {G. } ' <G . > 'dn (V . 16 ) 

0 1 3 



1 

/ {G. }<G .> 'dn (V.17) 

0 1 3 
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(V . 18 ) 



/ {G . } <G . >dn 

0 1 3 



/ {G. }<G. > ' {4 , }<G.>dn 

n 1 J J 



(V. 19 ) 



/ (G. }d n 

0 1 



(V. 20) 



To formulate these operators, the global shape function, G^, 
was defined on the local level by 



r - „ (i_1) m _ ( 1 ) 

G i - 9l © ^2 



(V. 21) 



where and g^ were defined by 



_ (e) 
g l 



(1 - y ~) for C in element (e) 
e 



for 5 not in element (e) 



(e) 
g 2 



for £ in element (e) 



for £ not in element (e) 
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and is the length of the eth element. The © notation 
in expression V.21 means that is the union of ^ 

and ^2 ^ * T * ie local (element) shape functions have the 
following properties. 



(i) 




0 if j ^ m 



(ii) g^ e) ( n j ) 



6 . 
13 



1 if i = j 

0 if i ^ j 



Having defined the local shape functions, the elemental 
matrix operators contributing to the global matrix operators 
V.16 through V.20 are. 



1 

/ {G.} <G . > 'dn 
0 1 J 



I 





1 

/ (G.}<G.>'dn ^ 
0 i 3 T 




1 



1 



(V.16 ' ) 



(V . 17 ' ) 
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2 



1 



£ 



/ {G . } <G . >dn 
J 0 1 3 



e 



(V. 18 ' ) 



6 



1 



2 



1 

/ {G. }<G .>'{'F}<G ->dn 





^1 

3 




(V. 19 ' ) 



/ (G. }dn | 
n 1 



0 



1 




(V. 20 ' ) 



The derivations of these operators are shown in Appendix E. 

Since there is extensive coupling in the field equations 
and there are three degrees of freedom at each nodal point, 
the numbering scheme for the system matrices was important. 
To minimize the bandwidth of the matrices, the numbering 
scheme represented in Figure 8 was used. 

Figure 9 represents the matrix A(t) of expression V.12 
for any three successive nodal points. The distribution of 
the elemental matrices for the FEM operators is shown, as 
well as the possible locations for the coupling of the field 
equations over an element. In addition. Figure 9 shows the 
bandwidth that would be observed for any n-1 element solu- 
tion. For this scheme, the bandwidth is nine. If coupling 
is not present, the bandwidth is seven. Figure 9 reflects 
the extensive coupling that is present in the model. All the 
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matrix elemental locations are filled except those which 
are shaded in. 

B. IMPLEMENTATION OF BOUNDARY CONDITIONS 

Having formulated the system matrices for the field 
equations, treatment of the boundary conditions will now be 
discussed. Each field equation is considered individually. 
1 . Fiber Heat Transfer Equation 

The fiber heat transfer boundary conditions in non- 
dimensional form from Appendix A are 



n = 0 A 



1 ft ■ V 0 ’ 11 + T^g -*»> 

QO 



(V . 22 ) 



n = 1 



1 It ■ -V 6 - 1 * - T^ig-^ 



(V. 23) 



Since the first term in expression V.13 is 



{G i }A 1 <Gj> ' {0} 



or in analogous form 



A 



1 




(V . 24 ) 



natural boundary conditions V.22 and V.23 may be directly 
substituted in equation V.24. 
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The response dependent parameters, h^, h 2 , T g , 
changing with time, are evaluated at the previous time 
step. Thus, the boundary conditions are incorporated in the 
system matrices as follows. 



( 1 ) 



-h^L: added to the stiffness matrix A(t) 



at location ^ 

(2) h n L - X|^-(T^ - T^ ) : added to the excitation 

1 T Q 00 

oo 

vector F (t) at location F. 



(3) 



(4) 



added to the stiffness matrix A(t) at 
location A 3n _ 2j3n _ 2 

h«L - - T^ ) : added to the excitation 

^ x cr °° 

no 



vector F (t ) at location F 



3n-2 



2 . Internal Flow Heat Transfer Equation 

The non-dimensional boundary conditions presented in 
Appendix A for the air energy equation are 



n = 0 r = 1 



(V. 25 ) 



9T _ 80 

8 n 8 n 



(V . 26 ) 



The essential boundary condition at ti = 0, r(n = 0) = 1 is 
imposed in the Galerkin equation as follows . The A 0 . row 
of the A (t ) matrix, the B 0 . row and the B. 0 column of the 
B matrix, and the F 2 location of the excitation vector, F(t), 
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are all set equal to zero. The B 2 2 location of the B 
matrix is then set equal to one. The natural boundary con- 
dition at n = 1 for the air heat transfer equation was treated 
in the same manner as that of the fiber heat transfer equation. 
Foregoing the individual steps the boundary condition is 
implemented by 



( 1 ) 

( 2 ) 



h 2 Lv l 



added to the stiffness matrix A(t) 
at location A 



h 2 Lv l 



aev l "4 *4 

±.(T’ - T ) 

T V g 

00 I 



3n-l , 3n-2 

added to the excitation 



vector F (t ) at location F^ . 



3 . Oxygen Transport Equation 

For the oxygen diffusion equation, the boundary 
conditions in non-dimensional form from Appendix A are: 



n = 0 |1 = U2 (4 - p) (V. 27) 

n = 1 f4 = 0 (V . 28 ) 

a n 



Since these are both natural boundary conditions, they were 
substituted for the first term in expression V.15 at n = 0 
and n = 1/ respectively. Coefficients and u> 2 are evalu- 
ated at the previous time step. The boundary conditions are 
implemented by adding 
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( 1 ) 



to the stiffness matrix A(t) at location 




and 

(2) w-p: to excitation vector F(t) at location 

F 3 

This concludes the discussion for the implementation 
of the boundary conditions. A word of caution is in order. 
After each time step integration, the time -dependent coeffi- 
cients of the boundary conditions ( i.e. , h^, A, v, w) 

are reevaluated. Before incorporating the updated coefficients 
into the stiffness matrix A(t), the previous values must be 
subtracted out. The results of not taking this into account 
will be obvious. 

The implicit system of ordinary differential equa- 
tions was integrated numerically by a modified implicit-Gear 
method developed by Franke [2] . The time-dependent coeffi- 
cients of the FEM operators in expression V.13, V.14, and 
V.15 were updated at the previous time. The reaction rate 
term (III. 26) was also evaluated at the previous time, and 
appears in the excitation vector F(t). In this way, the 
final system of ordinary differential equations was 

Aft)'}' + F (t ) = B l 

where A and B are (3n x 3n) matrices of temperature dependent 
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coefficients. F is a (3n * 1) vector arising partly from the 
reaction rate terms and partly from the boundary conditions. 
As noted, the elements of F are dependent on both temperature 
and oxygen concentration. 
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VI. RESULTS AND CONCLUSION 



A. NUMERICAL CONSIDERATIONS 

The input parameters which were varied in the computer 
analysis include (1) plate thickness, L, (2) fiber diameter, 
d, and (3) air flow rate, U w . Other parameters which could 
be varied are the ply thickness, D, ambient temperature, , 
and fiber emissivity, e. One set of initial conditions was 
used for the analysis. These are shown in Figures 10 and 11. 
For this initial effort, the actual initial conditions were 
not of prime consideration. The selection of initial condi- 
tions of Figures 10 and 11 will be discussed in Section VI. C. 
The finite element program calculates the remaining system 
parameters such as permeability, porosity, pressure differ- 
ential, and the response variables as functions of time and 
position. Parameters which are functions of temperature, 
such as k^, kg, k^, h^, p^, u^, and y are continuously updated 
during the transient analysis. 

The preliminary solution effort showed the system of 
equations to be very stiff (refer to Shampine/Gordon [24] 
and Gear [25] for a discussion of "stiffness"). Changing 
the integration algorithm from a sixth-order Runge-Kutta 
( IMSL subroutine DVERK) to a modified implicit-Gear method 
for stiff systems, developed by Franke [2], resulted in a 
significant reduction in CPU time. The computational effort 
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was carried out on an IBM 360/70. Typical runs required 20- 
30 minutes CPU time for problems for which ignition occurred, 
and 6-25 minutes CPU time for extinction problems. Ignition 
problems were terminated when the temperature at any nodal 
point exceeded 2200 degrees Fahrenheit, or when the graphite 
at a nodal point was totally consumed. Extinction problems 
were carried out to steady state. 

A twenty five nodal point model (75 o. d. e.) provided 
results which differed less than five percent from the re- 
sults of a 32 nodal point model (96 o. d. e.) and was adopted 
for all computer runs. For the twenty five nodal point 
model, approximately 275K bytes of core was required. The 
computer program which includes the FEM formulation and the 
integration routine as well as a sample input file is pre- 
sented at the end of the appendices. 

The numerical solution produced satisfactory results with 
one exception which will now be discussed. For a particular 
range of temperature and oxygen concentration (approximately 
1500 degrees Fahrenheit and .001 lbm/ft^, respectively) , 
there is a significant increase in reaction rate. This 
accelerated reaction rate produced a negative/positive 
oscillation for the nodal values of oxygen concentration. 

The oscillation occurred in a region of the plate where the 
oxygen appeared to be totally consumed. As discussed by 
Frank-Kamenetskii [3] , nth order reactions may yield zero 
values for the oxygen concentration and gradient within the 
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medium. The cause of the oscillation was as follows. 

Although the reaction rate is a function of temperature and 
oxygen concentration, numerically, it is treated as constant 
during an integration interval. As a result, the concentra- 
tion in the region of high reaction rate becomes negative. 

The reaction rate was updated for the next time interval by 
setting negative concentrations to zero. During this time 
interval, for the region of zero concentration, there is no 
reaction. Without consumption, the oxygen concentration will 
increase and the cycle repeats itself. The oscillation was 
aggravated by large plate thickness and low permeability. 

The instability occurring for the high reaction rate may be 
corrected by (1) numerically integrating the rate term in 
the time domain, (2) iterating until convergence is obtained 
between consecutive values of concentration, or (3) decreasing 
the time step and the length of the elements. As a matter 
of expediency, measure (3) was used to minimize the insta- 
bility. The fiber and air temperatures were essentially 
unaffected by the oscillations in the oxygen concentration. 

B. RESULTS AND OBSERVATIONS 

Table 1 presents the results of fourteen problems. In 
Table 1, the parameters which depend upon temperature (u , 

ir 

Re., and h.) are evaluated at ambient conditions for com- 
i l 

parison purposes only. These parameters, in fact, will vary 
during the course of the transient analysis. In all cases. 
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the initial condition on graphite temperature was 1050 
degrees Fahrenheit. 

Figures 10, 11, and 12 show the transient behavior of 
Case 1. For the initial conditions shown, ignition occurred 
in approximately eight seconds. The reaction rate and oxygen 
concentration responses for Case 1 are typical of the problems 
for which ignition resulted. Figure 13 shows that the air 
and fiber temperatures for Case 1 do not differ significantly. 
This behavior is typical for most problems. However, for 
plates with high porosity (greater than .8) and subjected to 
relatively high flow rates, significant differences in the 
air and fiber temperatures do occur (see Figure 14 for Case 
9) . In a previous effort by Vatikiotis and Salinas [26] , 
linearly varying initial conditions were investigated. Although 
the reaction rate used in that investigation differed from 
the one used in the present investigation, the overall re- 
sults remain valid. The main observation was that ignition 
is less likely to occur for thin plates with the ambient air 
entering the hotter plate surface. Since temperature con- 
trols the reaction in the kinetic regime, cooler air entering 
the hotter plate surface enhances the heat loss. This reduces 
the fiber temperature, thereby decreasing the reaction rate. 

In effect, the cooler air enters the plate where it is most 
needed. 

In the present analysis, observations were made by vary- 
ine one input parameter ( i.e. , L or d or U m ) while keeping 
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all others fixed. The effects of each parameter on the 
behavior will be discussed individually. 

1 . Effects of Exterior Velocity 

As shown in Figures 15, 16 and 17 (cases 4, 1 and 
2, respectively) for an initial condition of 1050 degrees 
Fahrenheit, three distinct regions of combustion were ob- 
served (extinction-ignition-extinction) . These became appar- 
ent as U was increased from low velocities (10 knots) to 

oo 

higher velocities (120 knots). Vulis [4] discusses an experi- 
ment of a heated carbon rod with an air jet impinging upon 
it. For a certain range of flow velocities, ignition was 
observed. However, extinction did occur for velocities less 
than and greater than the velocity for which ignition occurred. 
The behavior of the carbon rod and that of the graphite mat 
is similar. This behavior may be explained by the Semenov 
model of Figure 1. For high U^, the internal heat transfer 
coefficient, h^, will be large as a result of an increase 
in the pore velocity. Since the slopes of the heat loss 
lines, q , increase with increasing h^ for a constant T^, 
the graphite temperature decreases from the critical ignition 
point I. This causes extinction for the higher exterior flow 
velocities, U^. In contrast, for low U , the slope of q will 

00 J6 

be small and ignition is more likely to occur. However, the 
heat generation curve, q^, changes position due to a decrease 

in oxygen concentration. This reduction in oxygen concentra- 
tion is caused by the lower induced flow rate through the 
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plate. The overall shape of is flattened such that the 
line crosses q^ behind the new critical point I. Thus 
extinction is observed for this region of low U^. 

2 . Effects of Fiber Diameter 

Varying the fiber diameter with the remaining input 
parameters fixed, also yields three regions of combustion 
(extinction-ignition-extinction) for a fixed initial tempera- 
ture. The results are shown in Figures 16, 18 and 19 (Cases 
1, 9 and 10, respectively) . The fiber diameter and ply 
thickness determine the permeability of the plate. Similar 
to U^, permeability affects the pore velocity ( i . e . , low 
permeability -*• low velocity) . In turn, this affects the 
convective heat transfer and the amount of oxygen entering 
the plate. These three regions of combustion are explained 
by the Semenov model for the same reasons as those of varying 

U . 

00 

Porosity is a measure of void space and denotes the 

space available for oxygen. Porosity is associated with the 

internal geometry of the medium, as is permeability. However 

they are basically different parameters. Permeability 

(hydraulic conductivity) is associated with convective air 

flow through the medium. Thus, an n-fold increase in the 

fiber diameter and ply thickness will not affect a change 

2 

m porosity, but will affect an n order change in the permea 
bility. For an n greater than one, the result is a decrease 
in pore velocity without a change in porosity. Porosity 
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determines the maximum oxygen concentration per unit volume 

of fibrous mat. This indirectly affects the reaction rate. 

Thus, for a given u , the lower the porosity, the less oxygen 

P 

there is for combustion. It appears that high porosity 
would enhance combustion by the presence of more oxygen. 
However, higher porosities also provide more fluid per unit 
volume resulting in greater heat transfer. This is observed 
in Case 9, Figure 18. 

3 . Effects of Plate Thickness, L 

For the given initial temperature of 1050 degrees 
Fahrenheit, varying the plate thickness yielded two regions 
of combustion (extinction-ignition) . Extinction was observed 
for thin plates (Case 2, Figure 17) . As shown by Figures 20 
and 21 (Cases 7 and 11, respectively), ignition was observed 
when plate thickness was increased. Plate thickness also 
affects the pore velocity. Decreasing L, increases the 
pressure gradient across the plate, thus increasing the pore 
velocity. The result is that extinction is more likely for 
thin plates because of the enhanced convection heat trans- 
fer resulting from the increase in pore velocity. Whereas 
pore velocity is inversely proportional to L, it is an exponen- 
tial function of U^. For the range of parameters in this 
investigation ( i.e. , .1" < L < 3.", and 10 knots < < 100 

knots) , an n-fold change in will produce a greater change 
in Up, then will an n-fold change in L. Compare Cases 1 and 
2 to Cases 7 and 8. 
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A characteristic feature of ignition for thick 
plates was observed. As the thickness increased, the igni- 
tion region (the spatial location) moved closer to the 
entrance surface of the ambient air. Compare Figures 20 
and 21. 

C. DISCUSSION 

At this point, discussion is made of ignition tempera- 
tures for porous graphite plates. The ignition tempera- 
ture for Case 1 was 1050 degrees Fahreheit. This tempera- 
ture was obtained by varying the initial temperature until 
ignition occurred. In all subsequent cases, 1050 degrees 
Fahrenheit was adopted as the initial condition, and we 
observed whether extinction or ignition occurred. The igni- 
tion temperatures given here are the result of taking n = -j 
in the nth order reaction rate. Significantly different 
ignition temperatures will result for other n. For example, 
for n = 1 (1st order reaction) , the ignition temperature for 
Case 1 was approximately 1600 degrees Fahrenheit. As Frank- 
Kamenetskii [3] has suggested values of n between 1/3 and 
2/3, the average value of 1/2 was adopted for this analysis. 
The ignition temperatures obtained using n = j agree favorably 
with the experiments of Fontenot [1] . For fiber combustion, 
ignition temperatures for porous graphite plates are on the 
order of 1100 degrees Fahrenheit. 

Consider the results of varying U ot . Extinction occurred 
for Case 2 and Case 4; thus, the ignition temperature 
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for each case is greater than 1050 degrees Fahrenheit. 

Since ignition does occur at 1050 degrees Fahrenheit for 
intermediate values of U (Case 1, Case 7), ignition tem- 

oo 

perature is not a monotonic function of U . The precise 
determination of ignition temperature for a specific value 
of air velocity can be obtained by a large number of com- 
puter runs. A similar discussion can be made for the igni- 
tion temperatures dependence on fiber diameter. The general 
shape of these functions is represented by the shaded region 
on Figure 22. These suppositions may only be valid in the 
limited range of parameters associated with the present 
investigation . 

It is interesting to observe that changing the initial 
condition of the oxygen concentration only affects the early 
part of the problem. The response of the concentration is 
fast compared to that of temperature. After approximately 
one second, the transient behavior is the same regardless 
of the initial condition on the oxygen concentration. 

A brief description of the Semenov model was given in 
Section II. A. That discussion was for a lumped parameter 
model. The actual combustion process is more complex since 
there are spatial variations in temperature and in oxygen 
concentration. The results of this analysis show that both 
kinetic and diffusion regimes may exist simultaneously in a 
porous medium. Hence, it is not reasonable to restrict an 
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analysis of combustion in a porous medium to the kinetic 
regime. 



D. CONCLUDING REMARKS 

The model has provided considerable insight as to the 
behavior of fibrous composites subject to combustion. 
Improvements may be realized by extending the model to 
include : 

(1) The combustion of the epoxy matrix ( i .e. , start 
problem with first stage of combustion process) . 

(2) The change in geometry of the porous medium due 
to combustion (affects p, m and z) . 

(3) A more complex reaction ( i ,e . , secondary combustion 
of gaseous by-products, such as carbon monoxide). 

(4) The effect of gaseous by-products on fluid properties. 

(5) Generalization to two and three dimensional models 
(results in anisotropic properties) . 

The results show that certain properties, such as plate 
thickness, permeability and porosity, have significant effects 
on combustion. Therefore, one could select these design 
parameters to improve the flame resistance of composite 
materials. This would be beneficial to the survivability 
of a composite structure. A comprehensive set of analyses 
could be undertaken to couple the combustion problem to the 
strength problem. In this way, one could achieve a design 
which maximizes strength and minimizes combustion. 
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Another use of the present model would be in the study 
of energy systems utilizing particulate fuels ( i.e. , coal, 
biomass) . In this case, the design parameters could be 
selected to maximize the performance of the system. 



75 



APPENDIX A 



NOND IMENS IONAL I Z ATION OF FIELD EQUATIONS 



The field equations and boundary conditions are: 



3 T 



h. z 



3 T 



fx'VV-J # 1 -TT^t'VV +4 Hr(T g'*' = p g c gTt L (III ' 28) 



3 3T a 

3x (k a TF ) " mc a 3 x 



3T h.z 3 T 

— - + -i-(T -T ) = p c — rr~ 
p g a a a 3t 



(III. 37) 



- u pM - - (| )r(T g'*> - It 



(III. 46) 



at x = 0 



3T 

(k +k )-r- 3 - 

g r 3x 



h, (T-T) + ae (T*-T*) 
X g 00 g °° 



(IV. 1) 



T = T 
a ° 



(IV. 3 ) 



r i± 

B 3x 



U U - p4> ) 

p 00 



(IV. 5) 



at x = L 



3 T 



(k +k ) „ 
g r 3x 



-h (T -T ) - ae(T*-T*) 
Z g 00 g 00 



(IV. 2 ) 



3 T 
a 

3 x 



3 T 

_2 

3 x 



(IV. 4) 
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3 <£ 

3x 



0 



(IV. 6) 



The nondimensional variables are defined as: 



0 

r 

$ 

n 



T /T 

g °° 


nondimens ional 


T /T 

a °° 


non dimen s io na 1 




nondimens ional 


x/L 


nondimensional 



fiber temperature 
air temperature 
oxygen concentration 
distance 



The time variable, t, will not be nondimensionalized. 

Using the temperature of the fiber, T , to demonstrate 
the technique of transformation, we have. 



It follows that. 



T = T 0 

g oo 



3T 

2 

3x 



= T 



30 
« 3x 



2 

3 T 



3x 



_g . 



= T 



A 

3X 2 



Using the chain rule to transform x to n , gives 



3T 

2 

3x 



= T 



30 3 n 
® 3 n 3x 
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Since n = x/L, the partial derivatives become 



3T 

_2 = 
3 x 


T 

» 30 
L 3 n 




and 






3 2 T 

a - 

3x 2 


T 2 

«> 3 0 3 n 

L . 2 3x 
9 n 




Therefore , 






3 2 T 

a = 

3X 2 


T 2 

» 3 0 

~T 2 

L 3 n 




Similar transformations 


for T and <p 
d 


yield, 


3 T T . _ 

a _ °° 3 r 


a 4> 


^oo 3 $ 


3x L 3n 


3x 


FT Fn 


T - sr 


3 2 <J> 


^00 3 2$ 


3X 2 L 2 3n 


3x 2 ^ 


r 2 a 2 

L 3 n 



Substituting these relations into the field equations and 
removing the constants from the differential operators gives 



- TT^T <9 - r) + 4Hr( V + ) 



30 

P C rr 
g g3t 



T ” i_ (k i I) - mCa - T °° i I + V (0 . r) 
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P c tt: 
a a3t 
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^ 00 9 , p 3 $ > U p T °° 3 $ ^oo T oo *^ U p 3 T ,1, , . 

l 2 3 n B 3 n L 3 n L *3T a *3n $ ~ ( f )r(T g ,4>) 



3$ 
^ 3 1 



Upon rearrangement, the equations become, 



k [ ( V k r ) H ] -TI^T (0 ‘ r) + ^ r( V* } = p g c g L ^ 

— = »aV 2 M: 



3 u. 2 

c a °° r 



= L 



2 3_$ 
3 1 



Letting , 



h . zL 

A 1 = (k g +k r } A 2 = TT^ 



2 

pi 



. _ AHL 

A 3 T 



A, = p C 1 / 
4 g g 



v. = k 
1 a 



v 0 = me L 

Z 3, 



h . zL" 
i 

V 3 P 



v. = p c L" 
4 a a 



u)^ = B 



a). = u L 

2 p 



3 U 

“3 = 



(Jj . = 



1 1/ 

f 4> 



w 5 = L‘ 



the field equations may be written as, 



3_m ii) 

3n U 1 3 n 



- x 2 ( 0 -r) + x 3 r(T ,*) = 



A ii 

A 4 3 1 



(A . 1 ) 



3 , 3T, 

3n V 1 3n 



% r 

v 2 3? + v 3 (e - r) 



= V 



3 r 

4 3 1 



(A. 2) 
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f-r(w, |^-) “ 7n “ ^3 if 4 * " a, 4 r ( T a ' < f>) = 

3 ri J. 3 ri 2 3 n o3ri 4 g 


“5 3t (A>3) 



Noting that the coefficients X, v, u are response dependent 
variables, equations A.l, A. 2, and A. 3 are considered in 
their final form in Section III. 



Substituting the nondimens ional variables 
boundary conditions gives, 


into the 


X 1 If = h l L(0 - 1) + ' T !> 

oo J 


x = 0 


r = l 


o 

II 

X 


3$ /* N 

“13^ = U 2 (*“P) 

and 


x = 0 


X 1 If = -h 2 L(0-l) - - T 


x = L 


3 r _ 39 

3 n 3 n 


x = L 


|i = o 

3 r | 


A 

II 

X 



The dimensional form of the absolute temperature will be 
retained for convenience. 
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APPENDIX B 



POLYNOMIAL APPROXIMATIONS OF THERMAL PROPERTIES 



Relations giving the viscosity and thermal conductivities 
of air for varying temperature were required. A simple way 
to obtain values for these properties is to fit empirical 
data with 2nd order Lagrange polynomials . 

The general form of the 2nd order Lagrange polynomials 
for thermal conductivity and viscosity are 



<V T a2> <T a- T a3> . <V T a3> ( V T al> 

(T.,-T_,) al <T a2 -T a3 ) <T a2 -T al ) 



‘al a2 al a3' 



^a2 



(T a- T al> (T a- T a2> 
(T a3- T al> (T a3- T a2 



r k _ 
) a3 



(B . 1 ) 



l V T a2> <V T a3» 



. . ______ . <V T a3> <T a- T al> 

(T_,-T_,) (T_,-T_,) M 1 (T _-T ,) (T _-T , ) 



al a2 al a3 



'a2 a3 a2 al' 



<T a- T al> ( V T a2> 
lT a3- T al' IT a3- T a2' “ 3 



(B.2) 



Choosing three points from a range of temperatures that 
would be representative of those observed during the analysis, 
the corresponding values of the properties are: 
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(subscript) Temp 

(F) 

1 0 

2 700 

3 1500 



, Btu . , lbm > 

l lbm-hr-ft ; ^ft-hr ; 



0131 


.0394 


0284 


.0765 


0423 


.107 



Substituting these values into expressions B.l and B.2 
gives, 

(T -700) (T -1500) (T -1500) (T -0) 

k a = (0-700) (0-1500) ( • 0131) + (700-1500) (700-0) ( * 0284) 

(T -0) (T -700) 

(1500-0) (1500-700) 

(T a -700) (T -1500) (T -1500) (T -0) 

y = (0-700) (0-1500) ( * 0394) + (1500-0) (1500-700) (,0765) 

(T -0) (T -700) 

+ a a ( 1Q7} 

(1500-0) (1500-700P *- LU n 



These expressions reduce to. 



k 

a 



-3.232 x 10“ 9 T 2 + 2.412 x 10" 5 T + .0131 

a a 



Btu 

lbm-hr-ft 



u 



-1.041 x 10 _8 T 2 + 
a 



6.029 x 10 




+ 



.0394 



lbm 

ft-hr 
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Both expressions give property values which are within 
five percent of the data for temperatures to 2000 degrees 
Fahrenheit. 



83 



APPENDIX C 



TRANSFORMATION OF THE REACTION RATE TERM 



The expression for the oxidation rate of graphite fibers 
taken from Parker and Hottel [15] is 



r 

g 



9.55 x 10 6 

ml/ 2 *02 

gk 



,-44000 > 
exp ( — 7 ) 



R T 



gk 



gm 

2 

cm -sec 



(C.l) 



where R is the universal gas constant (1.958 cal/gmole-K) 
and P Q 2 is the partial pressure of oxygen (atm.). The 
partial pressure may be represented in terms of the tem- 
perature, T, and oxygen concentration, <}>, by the Ideal Gas 
law in the form given by Kanury [16] 

p o2 “ ( a~° r a * * (c - 2 > 

where M a / M 0 2 = *9 is the ratio of the molecular weights of 
air and oxygen; R^ (= 53.34 ft.-lbf/lbm R) is the gas 
constant for air. Substituting M /M^ ^ nto equation C.2 
and converting P^ to atmospheric units the expression 
becomes 



P 



02 



4.25 x 10 



-4 



R* T <P 



atm. 
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2 2 

Converting r to lbm/ft -hr from gm/cm -sec, equation C.l 

/ 

becomes 



r 



g 



7.04 X 10 10 „ 

372 *02 

gk 



,-44000, 
exp(-^- 7? — ) 



R T 



gk 



lbm 

ft 2 -hr 



Substituting in for and converting the absolute tem- 

peratures from Kelvin to Rankine, the rate expression 
becomes , 



r 

g 



4.014 x 10 



R a t 



a/2 



,-39883, 
exp ( — ) 



lbm 

ft 2 -hr 



The absolute temperature at the fiber surfaces is T . Thus, 

A A 

T in the numerator becomes T , and the expression becomes 



r 

g 



4.014 X 10 7 R T 1/2 
A g 



♦ exp ( 



-39883 ) 

T 

g 



To transform the reaction rate to lbm/ft 2 -hr, y, the ratio 
of a fibers surface area to its volume is placed in the 
numerator of r . The expression for y is, 



y = 



1 

2 




2 

d 



(C . 3) 



where d is the fiber diameter, and the j factor accounts for 
the uncertainty of the actual amount of fiber surface after 
the epoxy has burned away. The reaction rate of the graphite 
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fibers per unit volume is expressed as 



r 




-39883 




lbm 



g 



To obtain the heat liberated per lbm of fiber, r is 
multiplied by AH, the enthalpy of formation of carbon dioxide 
per unit mass of graphite consumed. This value is approxi- 
mately 14,085 BTU per pound-mass of fiber and was obtained 
from Kanury [16] 



To obtain the oxygen consumption rate, r is multiplied 
by the inverse of the fuel-air stoichiometric ratio. 



Heat generation rate = AH r(T^,<}>) 
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APPENDIX D 



JUSTIFICATION OF THE DANKV7ERTS ' BOUNDARY CONDITIONS 



The Dankwerts ' boundary conditions for the oxygen diffu- 
sion equation are: 



d£ 

dx 




U " P+J 



x = 0 



d<j> 

dx 



0 



x = L 



Bischoff [23] presents a discussion of these equations. A 
brief summary is provided for completeness. 

Figure 10 shows the one-dimensional oxygen transport 
region considered in this analysis. 



Flow 



x = 0 



x = L 



Entrance 

Section 



V u i 



Porous 

Plate 



4>/ B, u 



Exit 

Section 



<J) 2 / B 2 ' U 2 



The concentration of oxygen, <j>, for each region is distin- 
guished by subscripts; B and u are the diffusivity and the 
velocity, respectively, for each region. The oxygen diffu- 
sion equations for each section are as follows. 
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0 



x < 0 



(D.l) 



B 1 d 2 (J) 1 

u, , 2 
1 dx 


- = 0 X <_ 0 (D.l) 


B d 2 <j> 
U p dx^ 


- r (T, <j> ) =0 0 ^ x _< L (D . 2 ) 


B 2 d 2 <j> 
u 2 d? ’ 


d<)> ~ 

=0 x > L (D . 3) 

dx — 



with the general boundary conditions. 



$ ( -co ) 


= ^ (D.4) 


1^ 

O 

t — 1 

O. 


= <l>(0 + ) (D. 5) 


P^(O-) 


B 1 d<J» 1 (0 ) B d(j> ( 0 + ) , 6 

p u dx - * (0 ' u dx <D - 6) 

1 P 


P4> (L - ) 


= <P 2 ( L + ) (D. 7) 


p<j>(L - ) - 


B d 4» (L~) _ ( + ) _ B 2 d< ^2 (L ) , g. 

P u dx ^2 L > u 0 SZ (D - 8) 

P 2 


<t>2 ( + “) 


= finite (D.9) 
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For convenience, the parameters u^, will be treated as 
constant. In addition, the porosity, p, has been included 
to account for the discontinuity that exist for the concen- 
tration. The discontinuity is caused by the oxygen from a 
totally gaseous environment to one that is partly occupied 
by solid. 

An analytical closed-form solution for this set of 
equations is not possible because of the nonlinearity of 
equation D.2. However, the solutions of equations D.l and 
D.3 are, 



<j>^ = A^ + A^ expCu^x/B^) (D.10) 

<}> 2 = A 3 + A 4 exp(u 2 x/B 2 ) (D.ll) 

Applying boundary conditions D.4 and D.9, the following re- 
sults are obtained, 

A. = <J> 

1 “ 



The solutions D.10 and D.ll become 

‘t’j = ^ + a 2 ex P( u i x / B i) (D . 12 ) 
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2 



(D . 13) 



It would be necessary to have the solution for the non- 
linear equation D.2 to solve for A 2 and A^ . However, the 
constants need not be known to continue with the analysis. 
From equation D.12, 



^1 < 0 ) * ^00 + A ' 



and. 



d^O ) 

dx 



A 

B 1 2 



Substituting these into equation D.9 gives 



®1 u i + 

P^ + pA 2 - p ^(g± A 2 ) = *(0 ) 



_B_ d<j> (0 ) 

u dx 
P 



Cancelling terms yields, 



P4> 



(0+) _ df(0. ) 

y u dx 



Rearranging, the Dankwerts 1 boundary condition at x = 0 is 



d<j> _ U P 



dx 



B 



U - p o 



x = 0 
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From equation D.13, and noting that is a constant, we 
have 

d<P 2 



and 



d (J > 2 (L + ) 
dx 



0 



Substituting these expressions and equation D.7 into equation 
D . 8 , 



_ g 

P*(L") - P ~ = P4>(L") - ^(0) 

P 2 

Upon cancelling terms, the second Dankwets ' boundary condi- 
tion becomes. 



dj> 

dx 



0 



x = L 



An important underlying consideration for using the 
Dankwerts * boundary conditions is that it simplifies the 
analysis since equation D.2 may be solved independently 
without haveing to consider entrance and exit regions. 
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APPENDIX E 



DERIVATION OF THE FEM OPERATORS 

In the section on the finite element formulation, the 
following five differential operators were identified. 



1 

/ {G ! } <G >dn (E.l) 

0 1 ^ 

1 

J { G . } < G >dn (E . 2) 

0 1 J 

1 

/ {G . } <G . >dn (E . 3 ) 

0 1 3 

1 

/ {G. }<G'. >{y}<G.>dn (E . 4 ) 

0 1 3 3 

1 

/ (G. }dn (E . 5) 

0 1 



where the G^ are the global basis functions . These opera- 
tors are constructed on the element level by introducing the 
corresponding element basis functions, g^. The global and 
element basis functions are related by. 
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G . = g 



(i-1) 



where g^ and < 3 ^ are defined by 




g 



(e) 

1 



g 



(e) 

2 



(1 - y~) for £ in element (e) 
e 

0 for 5 not in element (e) 

Y~ for £ in element (e) 

e 

0 for £ not in element (e) 



and is the length of the (e)th element. 

The derivation of the local elemental matrices according 
to the Galerkin method for the global operations E.l through 
E.5 proceeds as follows: 

For operator E.l, 



global 



local 




I 




<g i g 2 >d ^ 



Noting that. 



gi 



1 



l 



e 
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<3o = T- 



the elemental matrix becomes 



r<ir> -<in> 1 

e e 




► f 


.-<T-> 2 <f )2 - 

e e 




e 




- 



-1 



For operator E.2, 



global 



local 



/ {G i }<G!>d n | U [ / 



a i g, 
e y l 






<g i 



Substituting in the local shape functions gives , 



(i - f-y 

e 



e 



< (- 1 , (f ) >d ? 

e e 



Carrying out the operations gives 



-1 1 



1 



g^>ds3 



94 



(- — + JL) 
v £ £ ' 
e e 



- ( — =y) 



— -L) 

' z z ’ 

e e 



(^-o) 



2 £ 



2£‘ 



ac | | 



"I I"! 



-1 



For operator E.3, 



global 



local 



/ {G } <G . >dn 

0 1 3 



» ot / 



Z g 
e y 



1 

g 2 ) 



; g x g 2 >d ?] 



Substituting in for the local shape functions 



(i - f-y 

e 



<f> 

e 



<(1 - 



f) ( f )>d « 

e e 



the elemental matrix becomes 



r (1 - 2 ^ + f 2 > <#-- 72 > 1 


Q 


c JC c X, 

e e 


« > -f 


„ 2 , 2 




(— — — ) (1-2 — + — ) 

4 „ 2 iX £ „ 2 ' J 





1 

2 



1 

2 
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For operator E.4, 



global 



local 



/ {G . } <G ! > { y } <G . >dn fr 

0 1 3 W 



£ 

/ 

0 



g l 



g 2 



<g i g 2 > 



V 

f 



<g 1 g 2 >d s } 



Multiplying out the matrices, the expression becomes. 



( gi 9 l 9 i t i-l +9 l 9 l g 2 T i ) <9 l 9 i 9 2 T i-l + 9 l 9 2 9 2 T i> 



L(g 1 g{g 2 » i _ 1 + 9 1 g 2 ? 2 T i ) + 9 2 9 2 9 2 |, i > - 1 



a? 



Integrating each term separately. 



£ 

e 

/ g^g^s = 
o 




l 



3 



{ 9 l 9 l 9 2 d « 




1 

3 



J q g ig -9 2 ac = 




1 

6 
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$ Q g l g 2 g 2 d ^ 




1 

6 



/ g^g-Jg^dc 
o ± ± 



l 



6 



/ o gi g 2 g-d ? 



1 

6 



e e 2 , 

/ g-ig 2 g 2 ds — — t d£ _ ~ 

0 1 2 2 0 JI 2 J 

e 

*e £ e 

j g ig 2 g 2 d? = / 3 d£ = T 

0 -L 2 2 o 2 



Substituting the values into the expression 



<- I Vi + 1 V 



<- 1 ' , i-i + 1 v 



<- I f i-i + I V 



<- § T i-1 + I V- 



factoring out a - 1/3, the local matrix for the global 
operator becomes 
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1 , 

/ {G. } <G . > { f } <G . >dn fc - 4 

1 * ^ 



¥ • , - Y . 

(t i-l-V ^~~2 i) ' 



<^2 i > V - 



For operator E.5, 



global 



local 



1 A A e l g l 

/ {G. }dn fr U[ / 

0 1 0 



d5] 



Substituting in the local shape function, and integrating, 
the expression becomes 




This last operator is used for the excitation vector as 
described in the FEM formulation. 
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REACTION RATE 



TABLES AND FIGURES 




FIGURE 1 Semenov model of* combustion. 







99 






FIBERS 



FIGURE 3 Separating a differential volume of porous medium 
into respective volumes of fiber and air. 



FIBERS 



^cond 
^rad | 



AHr(Tg,<*>) 



^ q rad 



dx 
x+dx 



(1 -p)d V 

FIGURE ^ Energy balance on the fibers. 



'conv. 



AIR 



^cond 
x 
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’conv 




"^■ q c 0 nd 
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a 



x+dx 
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FIGURE 5 Energy balance on the internal flow. 
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'}'CV») 



FIGURE 6 Molecule-mass balance on the oxygen. 
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FIGURE 7 Linear shape functions for the Galerkin formulation 
of FEM. 



Element O) ( 2 ) 

Local nodal point J 2 _J ^ i , 

Global nodal point } 2 3 n-2 

'*!■&, %=Q, 

*rR %r 3 

FIGURE 8 Numbering scheme used in system matrices. 
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FIGURE 9 Shows coupling in system of p. d. e. 
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SUMMARY OF RESULTS * 
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For D=.005" and initial conditions of Figure 10 
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FIGURE 10 Fiber temperature response of Case 1 (U=50kt. s d=,004 

T O \ 



CASE 1 
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FIGURE 11 Oxygen concentration response of Case 1 (U=50kt. s d=.004 
L=.2 5"). 
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FIGURE 12 Observed reaction rate for Case 
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FIGURE 13 Flow temperature follows the fiber temperature. 
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FIGURE 14 Flow and fiber temperature difference at high porosity and 
flow rate (for Case 9» p=.87, U=90kt.)* 
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FIGURE 15 For Case 4, U=£0kt. ( d=.004", L=.25" (for effect of varying 
exterior velocity, U^, compare Figures 15, 16 and 17). 
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FIGURE 16 For Case 1, U=50kt. , d=.004" 
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FIGUHE 17 For Case 2, U=90kt. , d=.004" 
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FIGURE 18 For Case 9» U^50k"t* » d=.002", L=.25" (for effect of varying 
fiber diameter, d, compare Figures 16, 18 and 19). 
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FIGURE 19 For Case 10, U=50kt. , d=.0048 
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FIGURE 20 For Case 7t Uj^Okt. , d=.004", L=.5" (for effect of varying 
plate thickness, L, compare Figures 17, 20 and 21), 
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FIGURE 21 For Case 11, U=90kt. , d=.004", L=l. 
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FIGURE 22 Shaded region represents shape of ignition 
temperature curve as a function of U^or d. 
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COMPUTER PROGRAM "CARBO" AND SAMPLE INPUT FILE 
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